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In the present work, we highlight the significant efTect that the simplest beyond nearest neighbor 
interactions can have on two-dimensional dynamical lattices. To do so, we select as our case example 
the closest further neighbor, namely the diagonal one, and a prototypical nonlinear lattice, the dis- 
crete nonlinear Schrodinger equation. Varying solely the strength of this extra neighbor interaction, 
we see examples of (a) destabilization of states that start out as stable in the nearest neighbor limit; 
(b) stabilization of states that start out as unstable in that limit; (c) bifurcation of novel states that 
do not exist in the nearest neighbor case. These dramatic changes are first theoretically highlighted 
through an analysis of a reduction of the problem to a few excited sites and the associated set of 
fT^ I conditions that govern their existence and their dynamical stability. Then, they are corroborated 

numerically through fixed point computations, spectral analysis and nonlinear dynamical evolution 
' simulations. 
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I. INTRODUCTION 
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' The broad theme of intrinsic localized modes has attracted a significant volume of attention from a wide range 
^SJ . of communities over the past two decades Such interest stems from their emergence not only in the modeling 
and computation, but perhaps importantly in experiments of a wide array of themes. These include, but arc 
not limited to, arrays of nonlinear-optical waveguides 0], Bose-Einstein condensates (BECs) in periodic potentials 
Q, micromechanical cantilever arrays Josephson-junction ladders Q, granular crystals of beads interacting 
through Hertzian contacts @, layered antiferromagnetic crystals 01 , halide-bridged transition metal complexes 0], 
Ph , and dynamical models of the DNA double strand Q ■ 

Among the many models of discrete systems (i.e., nonlinear dynamical lattices) that have been proposed for 
these physical set ting s, one that has been a central point of attention is the so-called discrete nonlinear Schrodinger 
(DNLS) equation [lO|. Part of the intrigue in this model lies in its apparent simplicity (since it incorporates solely the 
prototypical characteristics of interest, namely nonlinearity and a discrete form of dispersion), yet its considerable 
^ ■ wealth of nonlinear wave solutions and phenomena. Furthermore, its relevance as a suitable approximation of optical 
(in waveguides) [l, [ll| and atomic systems (in optical lattices) and the particularly simple form of breathers due 
to its separability of space and time variables in the standing waves of the DNLS only add to its appeal. 

One of the many themes that have been considered in the context of the DNLS equation is that of long(er) range 
interactions. In particular, the interest in this theme concerns the effects of potential inclusion of interactions beyond 
' those of the purely nearest neighbors of the standard DNLS. In that context, numerous interesting possibilities have 
I been found to arise. For example, it has been shown that interaction strengths with sufficiently slow decay (over 
■ space) can give rise to bistability of fundamental solitary waves (centered on a single lattice site) [l^. This, in turn, 
• • ' may play a role in soliton switching [l^l ■ On the other hand, such interactions may be relevant for energy and charge 
^ I transport in biomolecules [l5| and polymers [T6j . as well as in waveguide arrays. In the latter setting, one possibility 



' for their relevance is near the so-called zero-dispersion point [17[. But, arguably, a more relevant example is that 
5_( . of zigzag-shaped waveguide arrays which have not only been theoretically proposed [3, but also experimentally 
implemented and demonstrated to be valuable for promoting localization instead of diffraction even in the linear 
regime [l9| . As an aside, it is relevant to note that quantum variants of the nonlocal DNLS equation have been 
studied by means of the Bethe ansatz [2^. Furthermore, the role of nonlocal interactions in BECs has also been 
proposed to be critical in stabilizing rather unusual spatially periodic states such as the 3-site period waveforms 
identified in a discrete, long-range DNLS-type model proposed for ^^Cr, in the presence of an optical lattice [2l| . 
Although, a general formulation can be provided [l^, H^l for the existence and stability of DNLS standing waves 
from the well-known anti-continuum limit (of vanishing coupling between adjacent sites [23|), the number of studies 
that tackle nonlinear waves in higher dimensional longer range settings is very limited, to our knowledge, and chiefiy 
restricted to fundamental solutions; see e.g. [l^. Naturally, with the computational and analytical tools that are 
presently available, it is relevant to seek a deeper understanding of the role of longer range interactions in these and 
other wave systems, and some recent studies have been aiming in that direction, such as the Klein-Gordon chain 
analysis of |25| . 

Such an improved understanding in a special and suitably tailored, but also as we will demonstrate quite rich in its 
phenomenology, case example is the scope of the present work. In particular, we will focus on the setting where solely 
nearest neighbor and next-nearest-neighbor (that is, diagonal in our two-dimensional lattice) interactions are present. 
In this realm, we will demonstrate that the beyond-nearest-neighbor interactions are a powerful controller of both 
the existence and also of the stability of solutions and consequently of the dynamical evolution of the system. More 
specifically, we will demonstrate that solutions (such as the discrete vortices of topological charge S = 1), which are 
robust enough that they can be observed in photonic crystal experiments (26| , can be destabilized even by arbitrarily 
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small beyond- nearest- neighbor interaction (of suitable sign). Moreover, we will show that other states which are 
unstable in the standard case (such as the vortices of topological charge S = 2; see for theory [l^, [13] and for recent 
experiments [2^) will be stabilized when the next-nearest neighbor effect is sufficiently strong. Moreover, we will 
identify special limits (such as the degenerate unit diagonal neighbor limit) whereby even solitary wave (non-vortex) 
solutions will change (or exchange) their stability. This gives rise to unusual bifurcation events, such as a double 
pitchfork scenario that we identify in what follows, as well as to the emergence of novel and previously uncovered 
branches that solely exist because of the strength of the beyond-nearest-neighbor term. 

We should add a note in connection to the experimental implementation of waveguide arrays. The typical scenario 
of relevance therein involves a next-nearest-neiglibor interaction which is weaker than the nearest-neighbor one. 
However, as Figs. 1 and 2 of Ref. [l^ show, it is certainly possible in Id zigzag settings to create a next-nearest- 
neighbor interaction which is stronger than the nearest-neighbor one. This can be done systematically by simply 
modifying the angle of the zigzag lattice. On the other hand, in 2d admittedly this possibility is harder to realize. 
However, as is discussed in [2^, |30|, in the latter setting as well, it is possible to tune the interactions in femtosecond 
laser written waveguides with elliptical shape, by tilting the elliptical waveguides. While we are not aware presently 
of a 2d setting where such tilting (or other technique) has been able to produce a dominant next-nearest-neighbor 
interaction, nevertheless that regime is certainly of theoretical interest and, given the rapid progress of corresponding 
experimental technology, may soon become experimentally relevant as well, hence it is also considered here. 

Our approach in what follows will be two-pronged and will be based on a theoretical analysis of this setting from 
the anti-continuum limit and an accompanying set of continuation/bifurcation and dynamical evolution numerical 
computations from that limit. This will enable us to obtain systematic information about the existence of solutions, 
and also about their linear stability and to test these predictions against the corresponding numerical computations. 
Finally, the results will also be corroborated with direct numerical simulations to illustrate the stable or unstable (as 
appropriate) evolution of our identified waveforms. After presenting the model and theoretical setup in section H, 
our analytical and numerical results and their comparison will be given in section HI. Finally, in section IV, we offer 
a summary of the distinguishing features induced by the beyond nearest neighbor interactions in this system and a 
number of associated conclusions and potential future directions for the extension of the present study. 



II. MODEL AND THEORETICAL SETUP 

The DNLS equation on the two-dimensional square lattice of interest herein has a standard form [To| , 
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f^C(f>m.n + \(t>m.n\ 4>m,n=0, (1) 



where e is the coupling constant, and the linear operator C typically assumes the form of the discrete Laplacian A2. 
The latter is defined according to A20„i_„ ~ 4'm+i,n + 4'm,n+i + 4'm,n-i + 4'm~i.n ~ ^4>m.n- However, for our more 
general case, we will assume that 

(2) 

(mi,ni)eNN (m2,n2)£NNN 

That is we will consider both the effect of nearest neighbors (whose set is denoted by NN) and that of the next 
nearest neighbors (denoted by NNN). Notice that in the above, motivated by the optical waveguides problem, we 
use the optical notation where the evolution variable (propagation distance) is denoted by z. Also, in the expression 
of Eq. (I2|), the onsite term oc 4>m,n has been suppressed. 

Looking for stationary solutions, in the customary form = exp(iAz)um_„, Eq. ([T} leads to the time- 

independent equation: 

^"^m.n ^" ^^^ra,n \'^m,n\ '^m,n — 0- (3) 

Without loss of generality, we can rescale A = 1. Furthermore, to somewhat simplify notation, we will also use the 
vector formulation involving 1 = (m, n). 

Our analysis will take advantage of the well-established anti-continuum limit [23| , in order to develop a perturbative 
analysis from there. In particular, in that limit (of uncoupled sites), the energy of the decoupled oscillators assumes 
the form: 

E,{u) = Y,\u,\'-\\u,\\ (4) 
1 

Now. introducing the coupling adds a term to the energy with 

Ei[u) = -^e^ Ju/ {u{ui, + uiul) (5) 
11' 
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where the 1/2 prefactor is intended to avoid double-counting. In the setting described above, the kernel of interaction 
is non- vanishing only for |1 — l'| = 1 (NN), in which case Ji.i' = 1 and for |1— l'| = \/2 (NNN), in which case Ji_i' = k. 

The general persistence conditions [lol . [3ll [3^ of solutions from the anti-continuum limit (which can be straightfor- 
wardly written as u\ = e'^' ) demand that the unperturbed wave corresponds to an extrcmum of the perturbed energy 
for the solution to persist. This necessary condition suggests that the gradient of Ei(u) evaluated at a solution with 
multiple excited sites u\ = e*^' should vanish. Direct calculation of this yields for each site 1 the solvability condition 

^ JiP sin(0i - 0iO = 0. (6) 
i'#i 

Moreover, in order for the solution to be stable, the corresponding extremum has to be a minimum of the effective 
energy Ei [l^ [s^l (see also H^])- More specifically, the eigenvalues 7^ of the Hessian of Ei (evaluated at the 
above solution of the conditions of Eq. (O), are in fact connected to eigenvalues A (bifurcating from zero, when e 
becomes non- vanishing) of the original lattice dynamical problem. The connection is given by A| = 2e7j, to leading 
order [IS,!!!,!!!. 

We will apply these considerations predominantly for the case of 4-site squares with ncxt-nearcst-neighbor inter- 
actions, where we can extract specific analytical conclusions from the corresponding algebraic persistence and the 
above stability conditions. Nevertheless, we will also briefly mention the interest in potential generalizations for the 
8 site contours containing e.g. ((-1,-1), (-1,0), (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1)). 

The theoretical results will also be corroborated by means of full numerical solutions. Exact solutions of Eq. 
([3]) will be obtained by means of a Newton method. Upon generating such stationary solutions, their stability is 
examined through spectral stability analysis. To this aim, a perturbed expression of the form 

(/>m,n = exp(iz)w,„,„ -I- 5 exp(iz) [a,„,„ exp(Ai) -I- exp(A*i)] , (7) 

is substituted into Eq. ([T]). Here u^.n is the unperturbed stationary solution, 5 is an infinitesimal amplitude of the 
perturbation; A denotes the corresponding eigenvalues (which arc real or complex in the case of instability). This 
leads to the following linear equation for the perturbation cigenmode, 

where M is the Jacobian matrix, 

/ dFk/du, dFk/du* \ 
- [ -dF*/du, ~dF*ldu* ) 

and Fi denotes the left hand side of Eq. (jSj and the string indices fc} = m + {I — l)n, I = 1,2,.., N, map 
the N X N lattice into a vector of length N"^. Numerical solutions were sought for with the Dirichlet boundary 
conditions at the domain boundaries. Notice that given the localized spatial nature of the considered solutions, 
we expect that our numerical observations, for the range of parameter values considered herein, are essentially 
insensitive to the precise selection of boundary conditions. To generate numerically exact stationary solutions, the 
fixed point algorithm was iterated until convergence (typically with a tolerance of 5 x 10~^). Upon convergence, the 
spectral analysis of the stationary solutions was performed. The results are typically shown for 21 x 21 site lattices. 
When the solutions are found to be spectrally unstable, direct numerical simulations are performed (typically with 
a fourth-order Rungc-Kutta scheme), in order to detect the dynamical evolution of the instability. 



(8) 



III. ANALYTICAL RESULTS, NUMERICAL RESULTS AND COMPARISON 

We start with a consideration of the 4-site square, arguably the simplest two-dimensional contour that encompasses 
in a fundamental manner the higher-dimensionality of our setting. In this case, and using the relative phase variables 
4>i ^ O2 — Oi, (p2 " 63 — 6*2 and 03 = 6*4 — 9^ (for our 4-site contour with phase angles 0i,...,4), we can derive the 
following algebraic persistence equations from Eq. (|6]) 

= sin(0i) + fcsin((/)i + (/)2) + sin(0i + 02 + 03) (9) 
= sin(02) + fcsin(02 + <?!'3)-sin(0i) (10) 
= sin(03) - fcsin(0i + 02) - sin(02)- (11) 

It is now possible to manipulate the corresponding equations to get the set of solutions available for the system. As 
an indication of how to approach this problem, we note that adding Eq. Q and Eq. (|lip. we obtain sin(0i)-|-sin(03) = 
sin(02) — sin(0i -1- 02 -1- 03), which upon subsequent use of double angle formulas results in the conditions: either 
sin(^i±^) = or cos(^i^) = - cos(^^^±|i±^). 

Analysis of the resulting trigonometric conditions yields the following branches of solutions. 
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1. The standard discrete vortex of topological charge 5* = 1. This is the solution with = (0, 7r/2, tt, 37r/2), and 
01 =02 = 03 = 7r/2. This solution is well-known [lol. [27t to be stable for the nearest neighbor model of fc = 
with a double eigenvalue pair (to leading order) Ai^..._4 = zt2ei, a double eigenvalue at (due to the phase or 
gauge invariance of the model) and a higher order eigenvalue [that was calculated in [23| as As^g = ±\/32e^''^*]- 
In addition to illustrating the persistence (at least to the considered leading order) of such a solution, we have 
computed the Hessian of the perturbation energy of Eq. (O and have obtained the theoretical predictions for 
the corresponding eigenvalues in the presence of the next-nearest neighbor interactions, parametrized by k. We 
find that Ai_...,4 = ±2-\/efcj is a double eigenvalue pair (to leading order), while the leading order prediction for 
the remaining four eigenvalues is 0. Two of these will stay at due to the above mentioned invariance, while 
the other pair will bifurcate to higher order. Yet, here we would still like to focus on the significance of the 
lower order nearest-neighbor effect. It is remarkable that for these structures (which are called super-symmetric 
in [l^ [l^l because the leading order - A oc y/e - docs not contribute to their eigenvalues), the NNN effect is, 
according to this prediction, the dominant one. In fact, we can use an arbitrarily weak (even infinitesimally 
small in comparison to NN) "negative coupling" to render the configuration unstable. It should be noted 
that in the spirit of diffraction management and diffraction engineering jssl . js^l , this is certainly a scenario of 
potential physical interest. Nevertheless, we will also examine below numerous cases where positive k may have 
interesting implications on nonlinear wave stability as well. 

A relevant example of the corresponding branch of solutions is shown in Fig. [T] Numerical computations have 
been performed for e = 0.001. The top left set of panels showcases the positive k scenario (of stability), while 
the top right ones the negative k scenario (of instability). In fact, we observe that the situation is even more 
complicated because apparently higher order effects lead this double eigenvalue pair to be complex (although 
its real part is captured almost perfectly by our analytical prediction). The clear destabilization of the latter 
case is illustrated further in the bottom plot of k ~ —1.5. The dynamical evolution suggests a symmetry 
breaking between the amplitudes of the 4 sites (which start out as equal) that subsequently leads to a nearly 
periodic exchange of power between the 4 principal sites participating in the vortical structure. 

2. Another interesting solution is the so-called out-of-phase configuration with 6 = (0, tt, 0, tt). In this case all the 
0's are equal to tt. This configuration is also well known to be stable close to the anti-continuum limit for the 
focusing nonlinearities considered herein [lol | . In the present case, its corresponding eigenvalues are found to be 
Ai_2 = ±\/8ei, A3_4.5,6 = ±2-\/e(l — k)i (a double pair), while finally there is a pair at the origin, as expected 
due to the relevant invariance. Notice that in this case an instability is predicted when the NNN interaction 
strength overcomes the NN one i.e., for fc > 1. We will return to this effect later in our exposition. 

For now, let us comment on the very good agreement of the above predictions with what is shown in the left 
panel of Fig. [2] for e = 0.005 (that is used hereafter). Furthermore, the case oi k = 1.5 is selected on the right 
panel to indicate that although the solution is stable in the NN limit, a sufficiently strong NNN interaction 
may destablize it, leading to an amplitude symmetry breaking exchange of power among the 4 principal sites. 

3. Another principal, yet highly unstable configuration of the square contour that is predicted from our solvability 
conditions to persist is that with 9 = (0,0,0,0). Here all the 0's are 0. We examine this case mostly for 
completeness (and also because of its potential relevance and stability for the defocusing case of e < 0). The 
corresponding eigenvalues here are given by Ai^2 — ±A/8e, A3_4.5^6 = ±2-\/e(l + k) (a double pair), as well as a 
pair of zero eigenvalues. It should be noted here that this result, once again found to be in excellent agreement 
with our numerical computations in the left panels of Fig. [31 suggests a partial restabilization of this branch 
when the double pair becomes imaginary for A: < — 1. On the other hand, the dynamics of the branch shown 
in the right panel of the figure, is interesting in its own right as it suggests a pairing of (0,0) and (1,0) in an 
oscillatory pattern and of (0, 1) and (1,1) in a similar pattern [although this appears to change for sufficiently 
long time scales] . 

4. The next example is that oi 9 ^ (0,0, 7r,0), for which 0i = 0, while 02 = 03 = tt. In this case, our analytical 
calculation of the eigenvalues yields (in addition to the null pair) a pair at Ai_2 ~ ±2y/ke, a separate one at 

A3_4 = ±\/2e\/k + y/8 + Wi and one at X^^ = ±\/2e\/—k+\W+W . A brief inspection of these eigenvalue 
pairs confirms that there should always be at least 1 real and positive eigenvalue associated with this solution 
(for the focusing case, under study). If fc < 0, then there is exactly one such eigenvalue, while if fc > 0, there 
are two real eigenvalues. Hence, it should always be unstable. This, as well as our detailed prediction for the 
dependence of the A's on fc, are very accurately refiected in the full numerical computations of Fig. The 
left panel shows a prototypical example of the state and its systematic continuation over fc, while the right 
panel illustrates the oscillatory (yet not clearly periodic) pattern of exchange of intensity, upon the amplitude 
symmetry breaking that signals the pattern's predicted instability for fc = 1.5. 

5. We now turn to the case oi 9 ^ (0,7r,7r,0), for which 0i = 03 = tt, while 02 — 0. Once again, the eigenvalues 
of the linearization can be computed to leading order yielding in this case a double pair at the origin, one 
pair which is real in the absence of beyond-nearest-neighbor interactions and becomes modified according 




FIG. 1: The top left quartet of panels shows the real and imaginary parts of the solution, and its spectral plane (Ar,Ai) 
with A = Ar + iXi of linearization eigenvalues -all for k = 1-, as well as the dependence of the key (theoretically predicted) 
eigenvalues as a function of the next nearest neighbor strength of k. In the subplot that shows the eigenvalue dependence on k, 
the blue solid lines represent the numerical results, while the red dashed ones the analytical prediction. Notice the very good 
agreement between the two. While the left panel is for fc > 0, the right panel shows the case of fc < (the specific solutions 
and spectral analysis shown are for k = —1). The bottom panel illustrates the evolution of the intensity of the 4 principal 
sites of such a vortex (for k — —1.5) over the propagation distance z. In all of these graphs also below, the blue solid line will 
correspond to the site (0,0), the red dashed line to (0,1), the green dash-dotted to (1,0), and the black dotted to (1,1) for this 
4-site contour. 




FIG. 2: The left top panels correspond to the case of the (0, vr, 0, tt) state with k = 1.5. The instability of this branch for this 
supercritical case is evident in the spectral plane of the linearization. The bottom left panel corresponds to the continuation 
over the NNN interaction strength k and the squared eigenvalue (A'^) zero-crossing corresponds to the destabilization of the 
branch. Notice that we will often use hereafter this diagnostic (A'^) in cases devoid of complex eigenvalues, as its zero crossings 
are characteristic of stability changes of the solution and indicative of potential bifurcation points. This destabilization is 
dynamically illustrated in the right panel for k = 1.5, again featuring for sufficiently long propagation distances an amplitude 
symmetry breaking and an intensity oscillation of the 4 principal sites. 
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FIG. 3: The left panel of the figure shows the continuation and stability analysis of the 9 = (0, 0, 0, 0) branch. The top panels 
show an example of this branch for k = 1.5, while the bottom panel shows the corresponding single and double eigenvalue 
pairs. The right panel shows the evolution of the four central sites of the configuration (in the same way as before) for k = 1.5, 
indicating their paired oscillations between (0,0) and (1,0) and separately (0,1) and (1,1) for a lengthy interval during the 
propagation. 
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FIG. 4: The left panel shows a case example of the 9 — (0,0, 7r,0) configuration and its corresponding spectral plane (top). 
In the bottom, the dependence of the associated eigenvalues (see text) on k is given. The right panel shows the evolution of 
the central site intensities for the unstable configuration at fc = 1.5. 



to X^fi — ±2-\/ e(l — fc), in the presence of the fc-dcpendcnt ncxt-nearest-neighbor interaction. Notably, this 
dependence leads to restabilization of the configuration for fc > 1. Finally, the 4th pair is imaginary for fc > 
(but can lead to -further- destabilization for fc < 0), according to Xt^s = ±2-ye(l -I- k)i. Fig. [5] confirms once 
again the excellent agreement of the theoretical predictions with the relevant eigenvalue results (see the top 
right panel) and the existence of the instability in the absence of or for sufficiently weak beyond nearest neighbor 
interactions; see top left and bottom left panels. On the other hand, it also confirms the dynamical stability 
for the case of fc = 1.5 in the top left and bottom right panels. In the latter the small perturbation leads to 
bounded oscillatory dynamics, instead of the unstable evolution of fc = (bottom left). 

6. It is especially interesting to note that the configurations 6 = (0,7r,7r,0) and 9 = (0,7r, 0,7r) in the limit of 
fc = 1 become equivalent to each other. This is a byproduct of the equal strength of interaction of each of the 
sites with any one of its 3 (nearest or next nearest) excited neighbors. In this special limit, for both of these 
configurations, each of the excited phases of "sees" a neighbor with the same phase and two neighbors with 
TT phase and each of the tt phase excited sites "sees" another tt and two phases, rendering the configurations 
equivalent. This is manifested also by the equality of their respective eigenvalues in the expressions given 
above (they share a triple pair of O's and one pair of \/8ei). Given the stability change of these configurations 
at fc = 1 [0 = (0,7r,7r, 0) transitions from instability to stability as fc increases through the unit value, while 
6 — (0, TT, 0, tt) transitions in the opposite direction], we expect a potential bifurcation of a new branch past this 
critical point. Indeed, this is precisely what the analytical formulas of Eqs. ([5|)- ([TT|) predict. More specifically, 
the longer range interactions considered herein are not only responsible for stability changes (or exchanges), 
but additionally lead to the formation of entirely new branches of solutions that would be absent in the nearest 
neighbor limit. Such a branch is given hy 9 = (0, cos~^(— 1/fc), 2 cos~^(— 1/fc), cos~^(— 1/fc)). In this case, 

(^1 = 02 = -03 = C0S~^(-l/fc). 
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FIG. 5: The top left panel shows a case example of the 9 — (0, 7r,7r,0) for the stable case of fc = 1.5 and the unstable one 
of = 0. The top right panel confirms this transition from instability to stability as k is increased through following the 
corresponding eigenvalues numerically (blue solid) and analytically (red dashed lines), in excellent agreement betwen the two. 
The dynamical evolution of the principal 4 sites is demonstrated in the bottom left panel for the unstable A: = scenario, 
giving rise to a periodic emergence of an asymmetric pattern in the intensity of the sites. Finally, the bottom right panel 
case of fc = 1.5 only leads to (small fluctuation amplitude) bounded oscillatory dynamics even when perturbed, confirming its 
predicted dynamical stability. 



Wc can use the Jacobian formulation to provide explicit analytical predictions of the corresponding eigenvalues 
in this case, as well. In particular, in addition to the standard pair of eigenvalues at the origin, the other 3 
pairs are non- vanishing; A3_4 = i-^Se/fci, while X^^ = ±2y/e{l — k'^)/ki and Ay^g = ±2y^e(l — k'^)/k. From 
the above, it is clear that among the two extra vanishing eigenvalue pairs (at /c = 1) of the configurations 
9 = (0, 77, TT, 0) and 9 = (0, tt, 0, tt), in this "double pitchfork" bifurcation, one always exits as real and one exits 
as imaginary. This is a degenerate pitchfork bifurcation because at the critical point, there exist two eigenvalue 
pairs at the A = 0, which for fc > 1 move in different directions. As a result, this novel configuration created 
solely by the beyond-nearest-neighbor interactions will generically be found to be unstable in its interval of 
existence. These theoretical predictions are fully confirmed in Fig. |6l In particular, very good agreement (with 
the above theory) is obtained for the 3 pairs of bifurcating eigenvalues in the left panel, and the instability of 
the configuration with fc = 3 (left panel) is confirmed in the direct numerical simulations of the right panel. 

7. Finally, we touch upon a branch of solutions that is predicted by the leading order expansion and can be 
obtained for the values of e used here (in particular, for e = 0.005), but which does not exist as an exact 
solution in the case of fc = and hence we believe does not exist here either. Nevertheless, we have not 
proved this rigorously, since the proof would necessitate resorting to sufficiently high order expansions. We 
only infer this from the need to lower our tolerance to obtain the relevant solution (and the rigoro us p roof of 
its non-existence in for the fc = case, despite its theoretical proposition in physical setups in |35l|). Such 
asymmetric vortices have the phase profile 9 = (0, a, tt, tt -I- a), with 0i = a, (/)2 = tt — a and 03 = a. At the level 
of the leading order reductions, it is predicted that such vortices have a double pair of zero eigenvalues, while 
the other two pairs are located at X^^ = zL2^/e^ycos{a) — k and Ay^s = ±2-ye-\/cos(a) + ki. It is noteworthy 
that if this branch was an exact solution, it would naturally generalize the vortex branch with a = tt/2 and 
the mixed phase branch 9 = (0, tt, tt, 0) for a = n. It is then also natural to expect the extra vanishing pair of 
eigenvalues of such branches, due to essentially the invariance of the branch with respect to a. Nevertheless, as 
was shown e.g. for a = n/2 and fc = 0, this extra pair does not stay put at the origin, but instead it bifurcates at 
a higher order. The relevant "approximate" branch of solutions is shown in Fig. [7]for a = it/8. We can see that 
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FIG. 6: The left panel shows the (unstable) configuration with phase distribution 9 = 

(0, cos^^ ( — l/fc), 2 cos^^( — 1/fc), cos^^( — 1/fc)) for fc = 3. It also shows the dependence of the numerical eigenvalues 
(blue solid lines) for this branch and their comparison to the theoretical predictions (red dashed lines). The right panel 
confirms the instability of the configuration of the left panel for fc = 3, given the observed (strong) amplitude symmetry 
breaking in the intensities of the four sites from their (unstable) equilibrium values. 




FIG. 7: A case example of the approximate asymmetric solution branch with a = vr/S, shown in the left panel for k — 0.5. Its 
eigenvalues appear to be in very good agreement with the theoretically predicted ones (as a function of k), yet the solution 
is only approximate given that we have had to reduce the fixed point iteration tolerance to "converge" to it. The right panel 
shows the dynamical evolution of the solution at the left showcasing its lack of dynamical robustness due to a nearly periodic 
exchange of power between each of the two pairs of excited sites. 



it is always unstable for the range of considered /c's. This instability does appear to lead to symmetry-breaking 
pairwise oscillations/mass exchanges between the excited sites in the right panel of the figure. Although it is 
predicted that the instability should disappear for high enough fc's (beyond k = cos(a)), we were unable to 
converge to the solution up to that value of k (for a = tt/S), even for the reduced tolerance of 10~^ for the error 
in the convergence to the solution of our fixed point iteration used for this branch. While this is not conclusive 
in any way, it may be suggestive since these problems did not arise for the cases ofa = 7r/2ora = 7r examined 
above. 

Finally, to illustrate the powerful nature of the beyond nearest neighbor interactions as a "controller" of not only 
the existence but also the stability of complex nonlinear wave configurations, in Fig. [51 we present a select example 
from a configuration of a larger (8-site) contour in the case of a vortex of topological charge S = 2. In the case 
of fc = 0, the stability of this vortex has been analyzed in (see also the earlier numerical investigations of [36j ) 
and it was found that it was always unstable due to a higher order eigenvalue (proportional to e). In the present 
case, including the beyond nearest neighbor terms, leads to a dominant order prediction of 5 eigenvalue pairs at 
the origin, a double pair given by Ai^2,3.4 = ±2^/eki and a single pair of As^e = ±^/8eki. In Fig. [51 we observe 
that for A: < 0, these dominant eigenvalues give rise to a strong destabilization (with oc e, i.e., stronger than 
the fc = case) of the coherent structure. The associated instability is also evidenced dynamically in the figure for 
k = —1.5. On the other hand, for fc > 0, the relevant eigenvalues at ©(i/e) are imaginary and we observe from the 
figure that also all the higher order eigenvalues cross the stabilization threshold of A^ = and become imaginary 
(the last pairs cross for k w 0.7). As a result, the increase of the next-nearest-neighbor interaction is responsible for 
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FIG. 8: The top left panel presents the case example of the vortex of topological charge S — 2 for the unstable case of fc = —1.5. 
The relevant dominant eigenvalues shown also can be discerned in the blow-up of the top right panel that showcases more 
clearly the lower order eigenvalues which become stabilized due to the effect of k. Nevertheless, as shown in the evolution of 
the bottom left panel, the instability for k = —1.5 cannot be avoided and is manifested in the asymmetric evolution of the 
8 sites participating in the vortex (split into two pairs of four sites, namely the lower left ones and the top right ones). On 
the other hand, the bottom right panel shows that fc > has the same beneficial effect for the higher order eigenvalues (they 
become stabilized, and all have crossed = by the value of fc « 0.7). Yet in this case, also the leading order eigenvalues (well 
captured by our theory of the dashed red lines) are now imaginary and hence the vortex of 5 = 2 is completely dynamically 
stabilized. 



the complete stabilization of the vortex of topological charge S ~ 2. In that light, we conclude, that not only are 
such beyond-nearest-neighbor interactions potentially responsible for the destabilization of states that were stable in 
the nearest-neighbor-interaction limit (such as the vortex of 5 = 1). They are also potentially responsible for the 
stabilization of unstable states of that limit such as the vortex of 5 = 2. Finally, as we illustrated above, they arc 
also responsible for the emergence of novel states (such as the branch in item 6 above) and of unusual bifurcations 
(such as the double pitchfork that we obtained at the limit of A; = 1). 



IV. CONCLUSIONS AND FUTURE CHALLENGES 



In this paper, we have focused our interest on discrete two-dimensional dynamical lattices of the DNLS type, which 
arc a prototypical model for a variety of potential applications, including waveguide arrays in nonlinear optics and 
BECs in optical lattices (in the supcrfluid regime) in atomic physics. The principal theme of the study was the 
role of bcyond-ncarest-ncighbor interactions on the prototypical results that are known and understood for the more 
standard and more extensively studied case of the nearest-neighbor interaction [lol . [27| . 

What was found was that such additional interactions may play a critical role in shaping the associated dynamics. 
This is evident by their ability to destabilize stable configurations (such as the 4-site vortex of topological charge 
S ~ 1) and also their potential to stabilize unstable configurations (such as the 8-site vortex of topological charge 
S ~ 2), for sufficient strength and suitable signs of these beyond- nearest-neighbor effects. These two above mentioned 
examples are perhaps particularly notable because they belong to the category of the so-called super-symmetric 
states of [l^. For such states (for which the relative phase between adjacent excited sites in the contour is n/2), 
the contribution to the linearization Jacobian of the solvability conditions that we used to compute the full problem 
eigenvalues, remarkably, vanishes. Hence, in the nearest neighbor limit the dominant eigenvalues of such super- 
symmetric states arise with A of 0(e) or higher. In that light, the inclusion of beyond nearest neighbor interactions 



yields an effect which is dominant to leading order (with cx fee) even when k is small. It is thus natural to expect 
that especially in such super-symmetric cases, these longer range interactions play the role of a powerful controller 
affecting the potential stability of the ensuing nonlinear wave states. 

In addition to that possibility, we obtained a series of "solitonic" solutions (without a vortex structure), which in 
their own right had some interesting stability modifications, for sufficiently strong beyond nearest neighbor interac- 
tions. There, too, we saw solutions like the (0,7r,0,7r) start out as stable for small k but become unstable for large 
values of k, and vice versa solutions like (0,7r,7r,0) which start as unstable but are stabilized as k increases. In fact, 
these two solutions participate in an intriguing double pitchfork bifurcation at the degenerate limit of equal nearest 
and next nearest neighbor interactions of fc = 1. Off of this limit, we observe the potential of the next- nearest- 
neighbor terms to produce novel states such as those produced in item 6 above. Hence, this powerful controller of 
the beyond nearest neighbor interactions is responsible for the emergence of previously unfeasible wave states. 

We believe that through this prototypical example, we have made the case for the substantial relevance and 
interest within the consideration of interactions that go beyond the nearest neighbor effects. This is especially so in 
particular configurations (such as the super-symmetric ones) where the nearest neighbor eifects are not discernible 
and hence higher order interactions are dominant even as soon as they arise. It is thus an interesting direction to 
try to appreciate their effects more systematically, for different kernels, such as Gaussian, or exponentially decaying 
ones [37| . On the other hand, it would seem especially interesting to try to generalize relevant consideration to higher 
dimensions and to 3-dimensional solitons, vortices and vortex cubes [lOl. l38l|. to obtain a systematic view of beyond 
nearest neighbor interactions in such settings, as well. Such studies are deferred to future publications. 
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